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ABSTRACT 


A general finite element procedure for obtaining strain- 
energy release rates for crack growth in isotropic materials is 
presented. The procedure is applicable to two-dimensional finite 
element analyses and uses the virtual crack- closure method. The 
procedure was applied to non-singular 4-noded (linear), 8-noded 
(parabolic), and 12-noded (cubic) elements and to quarter -point 
and cubic singularity elements. Simple formulas for strain- 
energy release rates were obtained with this procedure for both 
non-singular and singularity elements. The formulas were 
evaluated by applying them to two mode 1 and two mixed mode 
problems. Comparisons with results from the literature for these 
problems showed that the formulas give accurate strain-energy 
release rates . 


INTRODUCTION 


Two-dimensional finite element analyses are widely used to 
obtain stress - intensity factors and strain-energy release rates 
for cracked isotropic and orthotropic domains. Several methods 
are available to extract the stress - intensity factor K and 
hence the strain- energy release rate G from finite element 
results [1-5]. For isotropic materials, stress - intensity factors 
have been used to predict fatigue crack growth and fracture. 
However, for composites, idealized as orthotropic or anisotropic 
materials, with cracks or delarainations , the strain- energy 



release rate G has been found to be more convenient [6,7] for 
these predictions . 

The s train - energy release rate for a particular crack length 
can be obtained with two finite element analyses using the same 
model but with two crack lengths differing by a small amount. 

This procedure is not preferred since the boundary value problem 
needs to be solved twice. A variety of methods [8-10] which 
utilize only a single finite element analysis are available to 
obtain s train- energy release rates. One such method, based on 
Irwin's virtual crack-closure method, uses the stresses ahead of 
the crack tip and the displacements behind the crack tip. 

Ryb icki and Kanninen [8] used this approach to obtain a simple 
formula for the strain- energy release rate for a cracked, 
isotropic domain modeled with four-noded quadrilateral, non- 
singular elements. 

In the recent years, several crack-tip singularity elements 
have been developed [2,3,4,11,12] . These elements produce the 
required singularity at the crack tip and give accurate stress 
distributions with fewer degrees of freedom than non-singular 
elements. For the special singularity elements with the near- 
field solution built in expl icitly [12], the K and G values 
are obtained as a part of the solution. For the more popular 
singularity elements, like the quarter-point (QP) and cub i c 
singularity elements, K and G are usually obtained using the 
crack opening displacement (COD) method [2-5] . This procedure is 
attractive for isotropic materials but has not been applied for 
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cracks along interfaces in bimaterial plates and is difficult to 
apply for cracks in orthotropic or anisotropic plates. 

The formula given by Rybicki and Kanninen [8] is attractive 
because the G values can be obtained easily from a single finite 
element analysis. Similar formulas would be very useful for 
higher order and singularity elements. The purpose of this paper 
is to present a uniform procedure for calculating strain- energy 
release rates that will apply to non-singular and singular 
elements of any order. The development of this procedure is 
based on Irwin's virtual c rack - c losure method. The procedure 
will provide a set of G formulas that depend on the crack opening 
displacements and the nodal forces at and ahead of the crack tip. 
The procedure outlined here will also serve as a basis for the 
development of similar formulas in the analysis of cracked 
orthotropic or anisotropic composite materials. These materials, 
however, are not considered in this paper. 

First, Irwin's virtual crack-closure method is reviewed. 
Next, a general procedure for calculating strain-energy release 
rates that applies to linear and higher order elements is 
presented. The procedure is then extended to quarter - po int and 
cubic singularity elements. Formulas for mixed-mode 
configurations are also obtained. Next, simpler versions of the 
formulas for the singularity elements are presented. The 
modifications required when the crack faces are loaded with a 
uniform pressure distribution are also presented. The G formulas 
are evaluated by comparing the calculated strain- energy release 
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rates with reference solutions from the literature for two mode I 
and two mixed-mode problems. 


a 

Ai , A2 , A3 , A4 
b 


E 


El , E2 



, F 


SYMBOLS 
crack length 

constants in assumed Oy distribution 

half-width of a center cracked tension 
specimen or width of a singe edge notch 
specimen 

Young's modulus for a homogenous 
isotropic plate 

Young's moduli in the bimaterial 
isotropic plate 

nodal forces at the i^^ node in the x- 
and y-directions, respectively 


G , Gj , Gij 


total, mode I, and mode II s t r a in - ene r gy 
release rates 


h 

K 

[Q] 


r 


u , 


V 


u 


m » 


V 


m 


half-height of specimen 
stress -intensity factor 

transformation matrix conne c ting forces 

F or F and constants A. 

X . y . 1 

1 •'1 

distance behind the crack tip on the y*0 
line 

uniform remote stress or uniform crack- 
face pressure loading 

displacements in the x- 
directions , respectively 

relative displacements 

node m with respect to 

node i : u - u - u , ; 

m m i 


and y - 

u and V of 
the crack tip 

V V - V . 

m m 1 
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v(x) 


displacenent interpolation function of a 
finite element 


^R(r) 

X, y 

<^x • f <^xy 


relative crack opening displacement at a 
distance r behind the crack tip 

Cartesian coordinates 

Cartesian stresses 

Poisson's ratios for isotropic, 
homogeneous and bimaterial plates 

crack growth increment, equal to the 
length of the element nearest to the 
crack tip in the x-direction 


VIRTUAL CRACK-CLOSURE METHOD 

Figure 1(a) shows a crack tip in an infinite isotropic plate 
subjected to remote mode 1 type loading. The normal stress 
distribution ahead of the crack tip and on the y - 0 line is of 
the form 

Ai 

- -7« + A2 + . . . (1) 

If the crack extends from a to a + A , for infinitismal values of 
A , the crack opening displacements behind the new crack tip 
will be approximately the same as those behind the original crack 
tip. Then the work necessary to extend the crack from a to a 
+ A is the same as that necessary to close the crack tip from a 
+ A to a. Irwin computed this work as 


W - I f v(r) a (A-r) dr (2) 

^0 y 

where v(r) is the crack opening displacement at a distance r 
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behind the crack tip at a + A . He then obtained the strain- 
energy release rate as 


^ Lim W Lim 1 , v , * . , 

- A-O A - A-0 2A dr (3) 

Rybicki and Kanninen [8] used Eq . (3) in their finite 

element analysis. They considered models with 4-noded 

quadrilateral elements only. The strain- energy release rate 


equations from 

Ref. 

( 8 ] are 



S “ 

1 

2A 

(F (V - 

V^,)l 

(4) 

^ii 

1 

' 2A 

- 

1 




where F and F are the nodal forces at node i in the 

^i ^i 

X- and y - direc t ions , respectively, computed from elements I and 
J in Fig. 1(b). The terms uj^ and are the displacements at 

node k in the x- and y-directions, respectively. The formula s 
in Eqs. (4) are very attractive because the G values can be 
computed from a single finite element analysis. These equations 
were shown to give accurate G values. 

Simple formulas, like Eqs. (4), will be very useful for 
higher order or singularity elements. This paper develops a 
uniform procedure, based on Irwin's virtual crack- closure method, 
that will give simple formulas, for the strain- energy release 
rates for higher order and singularity elements. 
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PROCEDURE FOR CALCULATING STRAIN - ENERGY RELEASE RATES 


In this section, a general procedure for obtaining the 
strain* energy release rates from a single finite element analysis 
is presented. First, this procedure is outlined for 4-noded 
quadrilateral elements. The procedure is then extended to 8- 
noded parabolic and 12-noded cubic elements. This procedure is 
also used to obtain G formulas for quarter -point (QP) and cubic 
singularity elements. 

In this procedure, the following assumptions are made: 

1. The finite element idealization in the immediate 
vicinity of the crack tip is symmetric about the x - 0 
line and symmetric about the crack plane (y - 0) . 

2. The normal and shear stresses on the y - 0 line and 
ahead of the crack tip are assumed to have the 
classical square-root stress distribution as in 

Eq . (1). In the limit, as the element size in the 

analysis domain becomes smaller, the classical square- 
root distribution will develop near the crack tip. This 
is true regardless of the type of element used to model 
the problem. Therefore, the form of Eq . (1) is assumed 

a priori instead of using the stress distribution given 
by the particular finite element model. 

3. The functional form of the crack opening profiles, u(r) 
and v(r), on the y - 0 line is determined by the element 
shape functions. 
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The normal stress Oy distribution of Eq . (1) can be 

determined from the nodal forces at and ahead of the crack tip 
(on the y - 0 line). While nodal forces are available at several 
nodes on the y - 0 line, in this procedure only the nodal forces 
from elements nearest to and around the crack tip will be used to 
determine the Oy distribution. Then the normal stress Oy and the 
crack opening displacements v(r) are used in Irwin's virtual 
crack- c losure method to obtain the G formulas for various higher 
order non-singular and singularity elements. 

Non - S ingular Elements 
4-Noded Quadrilateral Elements 

Consider a finite element idealization with linear 
quadrilateral elements symmetric about the crack tip as shown in 
Fig. 2. Ahead of the crack tip (at node i) and along the y 0 
line, the stress Oy is assumed to have the form 


o 

y ” -T^ 

+ A2 

for 

X 

IV 

o 




(5a) 


- 0 


for 

X < 0 




(5b) 

where Ai 

and A 2 

are 

unknown constants . 

The nodal 

forces 

F 

^i 

F . F , 

F shown in 

Fig . 

2(b) can be 

thought 

o f 

as 


cons is tent 

nodal forces 

acting 

at nodes i 

. j . 1 . 

m 

due to 

a 

prescribed 

stress 

distribution 

of the form 

in Eq . 

(5) 



Conversely , 

using 

these 

nodal 

forces , one 

could calculate the 
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constants in the assumed stress distribution of £q . (5) ahead of 

the crack tip along the x-axis. 

The stress distribution in £q. (5a) is valid only In the 

immediate vicinity of the crack tip since it is obtained from the 
near field solution at the crack tip. Also, since there are two 
unknown constants in the assumed distribution of £q. (5a), two 

forces are sufficient to determine the constants. The two forces 
F and F , computed from elements I and J nearest to and 

around the crack tip, can be used. These elements are shown by 
solid lines in Fig 2(b) . 

The work done by the assumed stress distribution on the 
boundary displacements of elements I and J is equated to the 
work performed by the forces F and F on the 

Yi Yj 

displacements Vj^ and vj as 
\ <^y(x) v(x) dx - - ^ 

where v(x) is the displacement interpolation function of the 4- 
noded element and A is the element length. Since a tensile Oy 
stress causes forces in the negative y-direction, a negative sign 


appe 

ars in Eq . ( 6 ) . 

For 

this 

linear quadrilateral 

element , 

the 

v(x) 

on the ij si 

de is 






v(x) - [1 

■ 

^i ^ 

/X. 


(7) 

Subs 

t i tut ing Eq . ( 7 ) 

into 

Eq. 

(6) and equating the 

mul tipi i 

ers of 

Vi 

and Vj yields 
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A A 


1* ■ Jo ' 75 ^ *2l *1 ♦ < t )''; 


A A 


■ • Jo ‘75 * *2' -' * ( f ) Aj> 


Eq. (8) may be written as 


where 


{F) - [Q] (A) 


(F) - {F F 


(A) - {A]^Va A 2 A) 


[Q] - 


^ [ 


Substitution of the calculated values of A^ and A 2 from Eq . 
(9) in Eq. (5) expresses the stress distribution Oy ahead of 
the crack tip completely in terms of nodal forces F and F 

Now consider the displacements behind the crack tip. The 
relative crack opening displacement vr at any distance r from 
the crack tip can be determined by the nodal displacements of 
nodes i, k and k' and the element shape functions as 


Vg(r) - (1 - ?) vj + I (v^ . v^,) 


( 10 ) 


Now that the stress distribution ahead of the crack tip and 
the relative displacements behind the crack tip are known, 
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Irwin's approach can be used to calculate the strain-energy 
release rates as 


r - i- 

I A-»0 2A 



a (A 

y 


r) Vj^(r) dr 


(11a) 


Lin 1 
A-0 2A 

Integrating 


Jo 


£q • 



A 2 I [ (1 




(V, 


(lib) and using Eq. (9) gives 


Vk.)ldr 

(11b) 


Lira 
I “ A-»0 


2 ^ 

2A 



+ 




,)] 


( 12 ) 


Because the relative displacement at the crack tip vj^ is zero, 
Eq . (12) reduces to 


G 


I 


Lira 

A-»0 



1 


(13) 


A similar procedure can be followed to show that 

= II - “0 <“k - "k-)I> <!'•> 

The limits in Eqs, (13) and (14) suggest that the crack tip 
element needs to be small, and as smaller elements are used at 
the crack tip, the correct limits for Gj and 

approached. Numerical experimentation is required to determine 
the necesary mesh refinement and the limits in Eqs. (13) and (14) 
can be then dropped in actual computations. 
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Eqs* (13) and (14) are identical to those presented by 
Rybicki and Kanninen [ 8 ] and are widely used in the literature 
[13-15]. Equations (13) and (14) are also applicable if 
triangular elements (Fig. 2(c)), obtained by collapsing the 
quadrilateral elements, are used at the crack tip. Note that for 
this idealization, the forces F and F are the forces at 

nodes i and j computed from elements I, J, K and L 
surrounding the crack tip. 

8 -Noded Parabolic Element 

To obtain the G expressions for the 8 -noded element, the 
procedure outlined above is followed. For this element, because 
forces are available at three nodes, i, j, and k (see Fig. 
3(a)), the Oy distribution along the y - 0 line is assumed to be 

Ai 

% “ 7 ^ 


The parabolic shape functions for this element are 


v(x) 


(1 






+ (- - + 2-^) V 

A A k 


(16) 


As before, the forces at the nodes i, j , and k can be 
related to the constants , A 2 , and A 3 of the assumed Oy 
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distribution (Eq. (15)) as 

(F) - IQ) tA) (17) 


where 



(F) 


(F F„ 

F 








^k 





(A) 

- 

(a^Ta 

A 2 A 




and 


/ 









168 

35 

8 



[Q] - 

1 

216 


224 

140 

96 






28 

35 

36 



The s tra in - energy 

release rates are obtained by 

performing 

integrations 

as in 

Eq 

. (11a) 

and 

using Eq 

. (17) 

as 



1 

2A 



(V 


m 




+ F (v- 


.)] 



1 

2A 



(u 


u 

m 


> ^ ^ <"i 
j 


U^,)l 




Note that the forces at node k do not appear in these equations 
because the relative displacement at node i is zero. 

These equations are similar to those for the 4-noded 
element. Eqs. (18) are also applicable if triangular parabolic 
elements (Fig. 3(b)), obtained by collapsing the parabolic 
elements, are used at the crack tip. Equations (18) were used in 
references 6, 7 and 16 to obtain G values for delaminated 
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composite laminates. Recently, Krishnamurthy et. al [17], using 
a procedure similar to the one described in this paper, derived G 
formulas identical to Eqs , (18) for this element. 


12-Noded Cubic Element 

For this element, because four forces are available at nodes 
i, j, k, and 1 (see Fig. 3(c)), the Oy distribution along 
the y - 0 line is assumed to be 


o 

y 



+ A^/x + A^x 


(19) 


The cubic shape functions for this element are 


v(x) - [1 - li J .9 (?)'- (^(|)']v^ 

^ 19(^) - (^)(I)' 

+ 1^1 


( 20 ) 


2'^L‘ 


As before, the forces at the nodes i, j, k, and 1 can be 
related to the constants Aj^ , A 2 , A 3 , A 4 in the assumed Oy 
in Eq . (19) as 

(F) - [Q] (A) (21) 

where 

T 


(F) - {F, 


F.. ) 


(A) - (Aj^Ta A^A A^A^/^ A^A^) 


14 


and 


“ ■ TO 


544 

105 

32 

14 

720 

315 

144 

63 

288 

315 

288 

252 

128 

105 

96 

91 


The strain-energy release rates are obtained by performing 
the integrations as in Eq . (11a) and using Eq . (21) as 


G_ - ' -Trr [F (v -v ,) + F (v -v , ) + F (v -v ,)] 
I 2h y^" p p'^ y^' n n' ' yj^' m m* ' ^ 


and 


( 22 ) 


j h 


II 2a ‘‘x^'^P “p ' ' ' n n’' x,_ ' m m 


Note that the forces at node 1 do not appear in these equations 
because the relative displacement at node i is zero. 

The strain-energy release rate formulas for this element are 
very similar to those of the linear and parabolic elements. The 
formulas given in Eq . (22) are also valid if triangular elements 

(Fig. 3(d)), obtained by collapsing one side of the cubic 
elements, are used at the crack tip. 


Singularity Elements 
Quarter- Point Singularity Element 

Henshell and Shaw [3] and Barsoum [4] showed that square 
root singularities are produced at the crack tip if the midside 
nodes of an 8-noded element (Fig. 3(a)) are moved to the quarter- 
point positions as shown in Fig. 4(a). The quarter -point nodes 
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on the y - 0 line in this figure are nodes j and 1 . 

For this element, the Oy distribution ahead of the crack 
tip on the y- 0 line is assumed to be like that given in Eq . 
(15), The shape functions of this element involve square root 
terms [4] and are 


v(x) - {1 . 3 Tj +2 \ 

A A 


)v^ + {4 


A 


/ ^ 1 
4 - ) V 

A ' j 




2 


(23) 


Using this shape function, the forces at the nodes i, j, and 
k can be related to the constants , A 2 , and A 3 of the 

assumed Oy distribution of Eq . (15) as 


{ F) - [Q] {A) 


(24) 


where 


{F) and (A) 

are defined in Eq . 

(17) , 

and 



10 

0 

-1 

[Q] 

1 

1 

40 

20 

12 



10 

10 

9 


For non-singular elements the inversion of the matrix [Q] 
was not necessary to determine (A). However, for the singular 
elements this was not the case. The constants A 3 , A 2 , and A 3 
are related to the nodal forces through 


{A) - [Q]-’- (F) 


16 


(25) 



where 

18 -3 6 

-72 30 -48 

60 -30 60 

The strain- energy release rates were obtained by performing 
the integrations as in Eq . (11a) and using Eq . (25). 



1 






where 


, . 3 3ir _ 

'll • • 


- 52 

33^ 

^2 

^ 2 

t . ^ * 2i» . 


17 

2 Itt 

*^21 2 16 ’ 

t22 

4 

t fl 21»r . 

*'31 ® “ ' 

f- _ 

- 32 

+ 21i- 

32 



(26a) 


(26b) 


A similar equation was obtained for Gji, where Fy is replaced 
with Fjf and v is replaced with u. 

In contrast to the regular parabolic elements, the strain- 
energy release rate equations for quarter-point (QP) elements 
have cross terms involving the corner and quarter -point forces 
and the relative displacements at the corner and quarter-point 
nodes. As before, the formulas given in Eq . (26) are also valid 

if triangular quarter-point elements (Fig. 4(b)), obtained by 
collapsing one side of the QP elements, are used at the crack 
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tip. 


Cubic Singularity Element 

Pu et al. [11] showed that a 12-noded cubic element has a 
square root singularity if the two side nodes are moved from the 
1/3 and 2/3 positions to the 1/9 and 4/9 positions, respectively 
(see Fig. 4(c)). This is analogous to moving the midside nodes 
to the quarter points for a parabolic element. For this cubic 
singularity element, the same procedure as above was applied to 
obtain the s tra in- energy release rate equations. 

For this element, the Oy distribution ahead of the crack 
tip on the y - 0 line is assumed to be as given in Eq . (19). 
The shape functions of this element, as is the case for the 
quar te r - po int element, involve square root terms [ 11 ] and are 


v(x) 


[1 - 


11 


J\ * 


?(^) 


3/2 


i 


(9 7f 


45 X 27 ,x, 

2 A 2 


3/2 




+ I- " y? + 18 * 


2 'A ■ *' A 


1Z(X) 

2 ^A^ 


3/2 


]v 


+ ? (^ ) 

2 A 2^A'' 


3/2 


]v. 


(27) 


With these shape functions, the forces at the nodes i, j, k, 
and 1 can be related to the constants A^, A 2 , A 3 , and A 4 
of the assumed Oy distribution in Eq . (19) as 

(F) - [Q] (A) (28) 


where {F} and (A) were defined in Eq. (21), and the new 


18 



transformation matrix [Q] is 


105 14 7 6 

315 63 0 -18 

315 252 189 144 

105 91 84 78 

The constants (A) can be related to nodal forces as 

(A) - IQ]'1 (F) (29) 

where 

216 -34 16 -54 

, / -1620 570 -300 810 

[Q]*'--47 

3240 -1290 1020 -2430 

-1890 770 -770 1890 

tm 

The strain-energy release rate equations were again obtained 
by performing the integrations as in Eq. (11(a)) and using Eq . 

( 2 9 ) a s 


r - ^ ^ 

I “ 2A 54 

(F 

^i 


+ 




^m' 

) J 

+ 

F 

''21<''p-''p'> 

+ 

'22<''n-''n-> 

+ 

'23<''.- 


)) 


F 



'32<''n-'’n'> 


'33<''«- 


) ) 

+ 

F 

^1 

"Al<''p-''p'> 

+ 


+ 

'a3<''»- 


))) 


(30a) 

where 
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ft 

- - 

11187 + - 

7155 
2 ^ 



*^12 

- 

38556 - 

24543 

2 

n 



- 

- 53055 + 

33777 

2 

7T 


^21 

- 

11396 

3 

9575 

8 

7T 


ti2 

- 

- 12936 + 

33003 

8 

n 


^23 

- 

17988 - 

45837 

8 

n 


^^31 


8453 
- 3 

3595 

4 

7T 


^32 


9804 - 

12411 

4 

n 


^33 

- 

- 13587 + 

17289 

4 

n 


Hi 

- 

6948 - 

17685 

8 

7T 


^42 

- 

- 23976 + 

60993 

8 

n 


^^43 


33372 - 

84807 

8 

' n 

A s 

imilar 

equation was obtained for Gjj 

replaced 

with 

Fx 

and V 

is replaced with 


(30b) 


formulas are applicable to triangular cubic elements (Fig. A(d)) 
obtained by collapsing one side of the cubic singularity 
elements . 


Mixed Mode Formulation 

The G formulas presented earlier apply for pure mode I and 
for pure mode II problems. In mode I (or mode II) problems with 
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parabolic elements, the force F (or F ) , at node k from 

*k 

element J (see Figs. 3a and 4a) is exactly equal in magnitude 

but of the opposite (or same) sign to the force computed from the 

element J', the symmetric element about the x-axis. In mixed 

mode problems, the deformation is neither symmetric or 

antisymmetric about the x-axis (or the crack line). Therefore, 

even if the finite element model is symmetric about the crack 

line, as shown in Fig. 5, the forces F ^ and F ^ will not 

^k 

be equal in magnitude and opposite in sign to each other. (See 

the appendix for details.) Because the force F for the 

^k 

regular parabolic element and F for the regular cubic element 
(see Fig. 3) do not contribute to G, the formulas given for 
these elements in Eqs. (18) and (22) are also valid for mixed 
mode conditions. 

However, the G formulas for the singularity element (Eqs. 
(26) and (30)) need to be modified for mixed mode conditions. 

The modification is necessary because, for quarter-point 
singularity elements, the nodal forces at node k have products 
involving the displacements at nodes 1 and m (see Fig. 4(a) 
and Eq. (26)). Similarly, for the cubic singularity elements, 
the nodal forces at node 1 have products involving 
displacements at nodes m, n, and p (see Eq . (30) and Fig. 

4(c)). The required modifications are explained below. 

In mixed mode conditions, the upper and lower crack faces do 
not deform symmetrically or antisymmetrically with respect to the 
crack tip. To close each of the crack faces from a+A to a, 
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different amounts of work are needed. The amount of work needed 


to close each crack surface can be found by computing the 
relative displacements of nodes on the crack surfaces relative to 
the crack tip and using the procedure outlined earlier. 

The G formula for the quarter point element, then, is 




(31) 


where ^12 » * • * ^32 defined in Eq . (26b). v^ and v^ , 

are the relative displacements of nodes m and m', respectively, 


with reference to the crack tip node i: v -* v - v. and 

m m 1 

T B 

V * V »- V. . The forces F and F are the forces 

m' m 1 y^ 

computed from the elements J and J', respectively (see Fig. 

5) . 

A similar equation was obtained for Gjj where Fy is 

replaced with F^ and v is replaced with u. For pure mode I 

B T B T 

(or mode II) conditions, F - p and F - F 

^k ^k 

(or and — - F^ ), and Eq . (31) reduces to 

^k ^k ’^k ^k 

Eq. (26). 

For mixed mode conditions, the G formula for the cubic 
singularity element Is 


22 


1 1 
2A 54 





h2<''n-''n-> 

+ 


21<’'p-’'p-> * 

'22<''n-''n’> 

+ 


31<''p-''p‘> * 

'32<''n-''n-> 

+ 

'33<''.-''«->' 

il'^p * ^ 42 '^n 

* '^43’.' 


(32) 


* '■y ,"=4l''p' ♦ ‘^43''«''> 


where ^\ 2 > .... t 43 are defined in Eq . (30b). The 

T B 

forces F and F are the forces in the y-direction computed 

from elements J and J*, respectively (see Fig. 5). 

A similar equation for Gjj is obtained where Fy is 
replaced with F^ and v is replaced with u. For pure mode I 
or mode II cases, Eq . (32) reduces to Eq . (30). 


SIMPLIFIED G FORMULAS FOR SINGULARITY ELEMENTS 

The G formulas for i-h^ non-singular elements, Eqs . (13), 

(14), (18), and (22), are simple. In contrast, the formulas for 

the singularity elements are awkward, even for pure mode I or 

pure mode II conditions. The formulas are even more complicated 

for mixed mode conditions. The complexity of these formulas can 

be traced to the terms involving the forces at node k for the 

QP elements (Eqs. (26) and (31)), and at node 1 for the cubic 

singularity elements (Eqs. (30) and (32)). The G formulas can be 

simplified considerably if these forces can be approximated. 

For the QP element, three constants A^ , A 2 , and A 3 were 
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assumed in the Oy distribution (Eq. (15)) because the three 
forces at nodes i, j, and k (Fig. 4) were used. However, if 
only two constants and A 2 are retained in the Oy 

distribution, only the forces at nodes i and j are needed to 
evaluate these constants. This assumption requires the forces at 
node k to be dependent on the forces at nodes i and j , and 
will simplify the formulas considerably. Setting the constant 
A 3 to zero in Eq . (25) yields 


F 



(33) 


Substitution of this value of F into Eq. (26) or (31) gives 


where 
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2A 
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A similar equation was obtained for Gjj where Fy is replaced 
with Fx and v is replaced with u. Note that this equation 
is valid for pure mode I, mode 11, and mixed mode conditions. 

This equation is considerably simpler and easier to use than 
either Eq. (26) or Eq . (31). 

A similar procedure can be followed for the cubic 
singularity elements. In this case, only three constants are 
used in the expression for Oy . Setting constant A 4 to zero in 
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Eq. (29) yields 


11 11 

F -F-^F +-TT*’ 

j-i ” yj >'k 

Substitution of this expression for F 


yi 


( 35 ) 

into Eqs . (30) and (32) 


gives 
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A similar equation was obtained for Gjj where Fy is replaced 
with Fx and v is replaced with u. Again, this equation is 
valid for pure mode 1, mode II, and mixed mode conditions, and is 
considerably simpler and easier to use than either Eq . (30) or 
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Eq. (32). 

Eqs. (26), (30), (31), and (32) will be referred to as the 

consistent formulas. Eqs. (34) and (36) will be referred to as 
the simpler formulas. The accuracy of the consistent and simpler 
formulas will be studied later in this paper by applying the 
formulas to two mode I and two mixed mode problems for which 
reference solutions are available. 

MODIFICATION FOR PRESSURE LOADED CRACK FACES 

The procedure developed earlier is valid for calculating 
strain- energy release rates for self-similar crack growth when 
the crack faces are stress free. When the crack faces are 
pressure loaded, the procedure outlined above must be modified. 

To illustrate this modification, consider a central crack of 
length 2a in an infinite plate. 

Figure 6 shows the region very near the crack tip with a 
uniform pressure loading of magnitude S applied to crack faces. 
As before, the s t r a in - ene r gy release rate is related to the work 
required to close the crack from a + A to a. The closure of 
the crack faces from a + A to a can be divided into two 
parts. In the first part, the applied pressure between a + A 
and a is erased and, in the second part, the stress free crack 
faces between a + A and a are closed. In the finite element 
analysis, the first part is equivalent to addition of nodal 
forces consistent to the opposite of the applied pressure at all 
nodes between a + A and a. For the second part, the general 
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procedure described earlier is used. 

For a regular 8-noded element, the uniform pressure S 

-i 

acting between a + A and a translates into three consistent 

S 2 S S 

nodal forces of values ^A, — jA, and ^A acting at nodes 1, 

j, and k respectively (see Fig. 6). To erase this pressure, 
the opposites of these forces are applied to nodes i, j, and k 
as shown in Fig. 6(b), and added to the finite element computed 

forces F^ , F^ , and F^ at nodes i, j, and k, 

^i ^k 

respectively, (see Fig. 6(b)). Thus, the forces at nodes i, j, 
and k, to be used in Eq . (18), will be 



Eqs . (37) show that the correct forces to be used in the G 

calculations for a crack face loading are sum of the computed 
finite element forces and the negative of the consistent nodal 
forces for the pressure distribution. The consistent nodal 
forces for an uniform pressure distribution of magnitude S for 
the various elements considered in this paper are shown in Fig. 

7 . These forces are then added to the forces computed at the 
nodes. The corrected forces, F , F , etc., are then used in 

y i y < 
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the appropriate G equations presented earlier to calculate the 


strain-energy release rates. 

EVALUATION OF G FORMULAS 

The formulas presented in this paper are used to calculate 

s t r a in - ene r gy release rates for cracks in two mode I and two 

mixed mode problems. The results are compared to those from the 

literature. The percent error is defined here as 

Present result - Reference result 

Percent error « — — — ★ 100 

Reference result 

Unless otherwise specified, a Poisson's ratio of 0,3 was used for 
all configurations analyzed. 

In all the problems studied, the finite element idealization 
at the crack tip had triangular elements with straight sides. 
These elements were obtained by collapsing one side of the 
parabolic or cubic elements. These triangular elements are 
preferred over their rectangular or curved counterparts as 
suggested in Refs. 4, 11, 18, and 19. 

Mode I Problems 

Remote Loading 

Center -cracked tension specimen . - The first example is a 
center cracked tension (CCT) specimen with a c r ack - 1 ength - t o - 
width ratio of 0.8 (see Fig. 8(a)). From the symmetries in the 
problem, only one-quarter of the specimen was analyzed. Fig. 
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8(c) shows the finite element idealization for the 8-noded 
parabolic elements and the 12-noded cubic elements, ihe specimen 
was analyzed four ways: (1) regular parabolic elements 

everywhere, (2) regular parabolic elements with QP elements at 
the crack-tip, (3) cubic elements everywhere, and (4) cubic 
elements with cubic singularity elements at the crack tip. For 
each case, the strain- energy release rates were obtained using 
the formulas presented earlier. Table 1 presents the normalized 
forces and crack opening displacements (COD), and the strain- 
energy release rates. The four G values are within about 3 
percent of an accurate value from Ref, 20, The models with the 
singularity elements, as expected, yielded the most accurate 
results. For these cases, both of the simpler formulas yielded G 
values which are more accurate than their consistent 
counterparts . 

Single edge notched specimen . - The second example was a 
single edge notched specimen with a/b - 0.8 (Fig. 8(b)). Again 
the finite element idealization of Fig. 8(c) was used. Table 2 
presents the normalized forces, COD, and the strain- energy 
release rates obtained for all 4 idealizations. The strain- 
energy release rates obtained were generally about 6% lower than 
an accurate solution given in Ref. 20. Again the models using 
the singularity elements yielded more accurate G values than the 
models using the non-singular elements. As in the CCT case, the 
simpler formulas yielded G values which are more accurate than 
their consistent counterparts. 
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Crack Face Loading 

As pointed earlier, the forces used in the G formulas must 
be modified for loaded crack faces (see Eq . (37)). To verify 

this modification, a center cracked plate with crack-face 
pressure loading was analyzed. Table 3 presents the forces on 
the nodes at and ahead of the crack tip obtained from the finite 
element analysis for the QP and cubic singularity elements. 

This table also presents the forces used in the G calculations. 
Because the s t r e s s - in t ens i ty factors and s train- energy release 
rates are identical for crack-face loading and remote loading, 
the COD should be identical for these two loadings. Therefore, 
the modified forces in column 4 of Table 3 should be identical to 
those obtained with the remote loading. A comparison of col umn 4 
of Table 3 and column 2 of Table 1 shows that this is true. This 
verifies the modification outlined earlier for the crack face 
pressure loading. 


Mixed Mode Problems 

To evaluate the accuracy of the formulas for mixed mode 
problems, an angle crack problem and a cracked interface problem 
were analyzed. The calculated strain-energy release rates are 
compared to those from the literature. 

Angle Crack 

Fig. 9(a) shows the configuration used for the angle crack 
in a finite plate. Fig. 10(a) and (b) present the finite element 
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idealizations with 8-noded parabolic elements and 12-noded cubic 
elements, respect ively ^ Fig. 10(c) shows the details at the 
crack tip. With the parabolic elements, quarter point 
singularity elements with A * 0.04a were used at the crack tip. 
In the cubic element model, cubic singularity elements with A » 
0.04Sa were used at the crack tip. 

Fig. 11 shows the forces at the crack tip. The forces at 
nodes i, j, and k and the displacements at nodes m and 1 
are needed to compute G. These forces and displacements, 
represented as components that are normal and tangential to the 
crack plane, can be computed as follows. First, the Cartesian 
forces and Fy at nodes i, j, and k are obtained. At 

node i (Fig. 10), the forces F and F can be computed 

from elements I, J, K, and L. The forces F , F , F , and 

X . y, X,’ 

J k 

F can be computed from element L alone. The tangential and 


normal forces at 

the 

crack 

tip, F^ and F , 

t . n . 

1 1 

can be computed 

as 





1 

F 

X . 
1 

cos ^ 

-f F sin^ 

^i 

(38) 

F 

n . 
1 

- F 

X , 
1 

s in^ 

+ F cos^ 

^i 



with similar equations for the forces at nodes j and k. 
Similarly, the normal and tangential displacements at nodes m, 
m* , 1, and 1’ relative to the crack tip can be computed by 
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1 


first calculating the Cartesian relative displacements as 



(39) 

V - V - V . 

m mi 

A similar set of equations were used for nodes m' , 1, and 1' . 

The normal and tangential forces and displacements are then used 
to obtain the Gj and Gjj values. 

Table 4 presents the Gj and Gjx values obtained with 
quarter -point and cubic singularity elements. For comparison, 
this table also contains the results from Wilson's collocation 
procedure [21] and Tan et al . ' s boundary force method [22]. 

In view of the slight (2 percent) difference between the 
collocation and boundary force method results, the present 
results are compared with the average of these reference results. 
The present results are in good agreement with the reference 
results. The accuracy is better for mode I than for mode II. 
Again the differences between the simpler and consistent formulas 
are negligible. 

Cracked Bimaterial Plate 

The second mixed mode problem analyzed was a biraaterial 
plate with a central crack along the interface (Fig. 9(b)). 
Erdogan and Gupta [23] analyzed this problem for an infinite 
plate with uniform pressure applied to the crack faces. Their 
calculated G values will be used to evaluate the present 
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results. 


Because of symmetry, only one-half of the plate was 
idealized, as shown in Fig. 12(a). Fig. 12 shows the finite 
element model using 8-noded parabolic and QP singularity 
elements. Table 5 presents the strain-energy release rates 
obtained using Eqs. (31) and (34) for three sets of material 
combinations. The strain- energy release rates given in Ref. 23 
were divided by a factor 2 because Erdogan and Gupta's definition 
of G was different from the classical definition (see Eq. (37) 
on p. 1097 of Ref. 23). Excellent agreement is obtained between 
the finite element computations and those of reference 23. 

In this bimaterial case the power of the stress singularity 
at the crack tip is of the form -l/2±ioc, where oc depends on the 
relative properties of the two materials. This indicates that 
the stresses oscillate near the crack tip. The solution in Ref. 
23 includes the oscillatory part of the stresses. In the present 
finite element analysis, however, this oscillatory behavior is 
neglected. The accuracy of the present results show that 
contribution of the oscillatory part to the strain-energy release 
rate is negligible. 


CONCLUDING REMARKS 

A general procedure is presented for obtaining the strain- 
energy release rate G for cracks using the virtual crack- 
closure technique in finite element analyses. The procedure was 
applied to non-singular 4-noded, 8-noded (parabolic), and 12- 
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noded (cubic) elements and then to quar ter - po in t and cubic 
singularity elements. The procedure assumes that the finite 
element idealization is symmetric about the crack plane at the 
crack tip and symmetric about a line normal to the crack plane at 
the crack tip. The procedure uses forces at and ahead of the 
crack tip and displacements behind the crack tip. With this 
procedure, simple formulas were obtained for the non-singular 
elements . 

For the singularity elements, two types of formulas - 
consistent and simpler - were obtained. The consistent formulas 
use three nodal forces for the quar t e r - p o i n t singularity element, 
and four nodal forces for the cubic singularity elements. In 
contrast, the simpler formulas use two and three nodal forces for 
the quar ter - point and cubic singularity elements, respectively. 

A slight modification was necessary to use the formulas when 
the crack faces were subjected to pressure loading. The finite 
element forces need to be modified to account for the pressure 
loading. The modification is simple, accurate and is easy to 
implement in the actual analysis. 

Two mode I and two mixed mode problems were analyzed using 
the non-singular and singular elements mentioned above. The 
strain- energy release rates for these problems were obtained 
using the formulas derived in this paper. Comparisons with 
accurate or reference solutions from the literature showed that, 
for the configurations analyzed, the G formulas yielded accurate 
results. For all the configuration analyzed, the simpler 
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formulas for the singularity elements yielded strain- energy 
release rates which were more accurate than their consistent 
counterparts . 

With the singularity elements, fewer elements (and, hence, 
fewer degrees of freedom) are needed to obtain accurate 
solutions* The simpler strain -energy release rate formulas 
developed in this paper for singularity elements considerably 
simplified the G formulas without sacrificing accuracy. These 
formulas will enhance the utility of these elements in the 
analysis of cracked isotropic and bimaterial configurations. 

The analysis and procedure presented in this paper should 
provide the basis for future development of formulas for strain 
energy release rates for cracked orthotropic and anisotropic 
materials . 
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Appendix 

Inequality of Forces at Node k 

This appendix will show that the forces at node k (see 
Fig, 5(b)) computed from the elements J and J' will not be 
the same for mixed mode conditions . 

A mixed mode condition can be thought of as the 
superposition of mode I and mode II conditions. Fig. 13 shows 
the forces at node k computed from elements J and J' for 
mode I, mode II, and mixed mode conditions. For the mode I (mode 
II) case, symmetry (antisymmetry) conditions about the y - 0 line 
require that the x- and y-forces be symmetric (antisymmetric). 

Note that the forces for mode I and F^^ for mode II will 

X y 

not, in general, be equal to zero, (The requirement of a zero x- 
force at node k for mode I and a zero y-force at node k for 
mode II can still be satisfied because elements to the right of 
elements J and J' contribute to these forces.) 

Superposition of the forces at node k from the mode I and 
mode II conditions will yield the forces for the mixed mode case 
as shown in Fig. 13. The y-force at node k from element J is 

F^ + F^^ while that from element j" is - F^ + F^^. Similarly, 

y y y y 

the x-force at node k from element J is F^ + F^^ and that 

X X 

9 I II 

from element J is F - F . Thus, the forces at node k 

X X 


37 



computed from the top and bottom elements will not be equal. 
Similar arguments hold for node 1 in Fig. 5(c). 
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Table 1 • Comparison of strain-energy release rates for a 


center 

cracked tension 

specimen subjected to 

remote 

uniform 



stress S 




(a/b- 0.8 ; A/a 

* 0.0625 ; plane strain ) 



G [20] - 3.298 S_ ?_ « ( 1 

Reference E 

) 


Element 

Nodal Forces 

Relative Crack -y- 

EG 



Percent 

Used 

Fy/Sa 

opening S 

wa(l-i^ ) 

Error 



Displacement 





Ev/Sa 



Non - S ingular 

Elements 




8 - noded 

i -0.3218 

m 2.179 



Parabolic 



3 . 197 

-3.0 

E lament 

j -0.2969 

1 1.487 



(Fig. 3(a)) 






1 2 - no de d 

i 

-0 . 2552 

P 

2 . 240 



Cubic Element 

j 

-0.2360 

n 

1 .835 

3 .239 

-1.8 

(Fig. 3(c)) 

k 

-0 . 1203 

m 

1 .269 



Sineularitv Element 






QP S ingular i ty 

i 

-0 . 1068 

m 

2 . 276 

3 . 242 

-1.7 

Element 

j 

-0.4348 

1 

1.156 

(3 . 269)* 

(-0.9) 

(Fig. 4(a)) 

k 

-0.0941 






Cub i c 

i 

-0.0806 

P 

2 .292 

3 . 271 

-0.8 

Singularity 

j 

-0 . 2406 

n 

1 . 544 

(3.285) 

('0.4) 

Element 

k 

-0 . 2471 

m 

0 .7767 


(Fig 4(c)) 

- 

1 

-0 . 0698 






Values in parenthesis were obtained with simpler formulas. 
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Table 2 - Comparison of strain-energy release rates for a 
single edge notched specimen subjected to remote uniform 


stress S 

(a/b - 0.8; A/a - 0.0625; plane strain) 

2 2 

S IT a(l-i> ) 

G [20] - 143.8 

Reference E 

EG 

Element Nodal Forces Relative Crack -= ***o Percent 

Used Fy/Sa Opening S^ira ( 1 ) Error 

Displacement 

Ev/Sa 


Non-Singular Element 


8 - noded 
Parabolic 

i 

-1.985 

m 

15.95 

128.9 

1 

O 

O 

Element 
(Fig. 3) 

j 

-1.414 

1 

10.19 



12 -noded 

i 

- 1 . 602 

P 

16 . 74 



Cubic 

j 

-1 . 346 

n 

13 . 10 

134 . 6 

-6.4 

Element 
(Fig. 3) 

k 

-0.4236 

m 

8 . 645 




S ineular i tv 

Elements 




QP 

i 

-0.7192 

m 17.08 

134.8 

-6.5 

S ingular i ty 




(137.4)* 

(-4.5) 

Element 

j 

-2.432 

1 7.804 



(Fig. 4) 







k 

-0.2309 





Cubic 

i 

-0.5052 

P 

17 . 37 

138 . 6 

-3.7 

Singularity 

j 

-1.574 

n 

10 .85 

(139.0) 

(-3.3) 

Element 

k 

-1.163 

in 

5 .191 



(Fig. 4) 

1 

-0.2596 






it 

Values in parenthesis were obtained with simpler formulas 
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Table 3 


- Nodal forces for a center cracked tension specimen 
with uniform crack face loading of magnitude S 

(a/b - 0.8; A/a - 0.0625; plane strain) 


Fy/Sa 


Forces 

computed in the 

Forces needed 

Forces to 

f ini te 

element analysis 

to erase the 

be used in 



crack face 

the formu- 



loading 

la* 



i 

-0.1068 

0 

-0 . 1068 

QP 

Element 


-0.3931 

-0 . 0417 

-0.4348 

(Fig. 4) 


k 

-0 . 0733 

-0 . 0208 

-0 . 0941 


Cubic 

i 

-0 . 0785 

-0.0021 

-0 . 0806 

Singularity 

Element 

J 

-0.2312 

-0 . 0094 

-0 . 2406 

(Fig. 4) 

k 

-0 . 2096 

-0.0375 

-0 . 2471 


1 

-0 . 0562 

-0 . 0135 

-0.0698 


Forces in this column are identical to those computed with a 
remote uniform stress of magnitude S applied on the line 
y - h (see Fig. 8). 
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Table 4 - Comparison of strain- energy release rates 
for an angle crack in a finite plate 
(a/b » 0.5; ^ - 45"; plane strain) 


Element 

Used 


EG. 


2 2 
S fta(l-i)^) 


EG 


II 


S^»a(l--j)^) 


QP 

Singularity 

Element 


1.A36 0.3378 

(1.447)* (0.3399) 


Cubic 

S ingular i ty 
Element 


1.433 0.3446 

(1.444) (0.3354) 


Reference Results 

Wilson [21]"*" 1.440 0.325 

(Collocation) 

Tan et al . [22] 1.416 0.3295 

(Boundary Force 

Method) 

Values in parenthesis were obtained with simpler formulas. 
“^These values were taken from Figs. 60 and 61 of ref. 21 
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Table 5 - Comparison of strain-energy release rates for center 
crack along the interface in bimaterial plate with 
crack face pressure loading S 
(a/b - 0.1; plane strain) 


El G 





S2 a ( 1 

- 7)2i) 



Material 

Mode I 

Mode II 

Total 

Total 

Percent 






Ref [23] 

^ Error 

1 

2 





Aluminum 

Epoxy 

10.35 

0.7262 

11.08 

10 .92 

1.42 



(10.43) 

(0.6296) 

(11 . 05)''' 


1.21 

Steel 

Epoxy 

29.56 

2 . 808 

32.05 

31.66 

1.23 



(29.77) 

(2 . 168) 

(31 . 94) 


0 . 88 

Steel 

Aluminum 

1 .934 

-0.0567 

1.877 

1 . 980 

-5.19 



(1 . 947) 

( -0 . 0567) 

(1.890) 


-4.56 


MATERIAL PROPERTIES 


Material 

E. 

ps i 

■i> 

Epoxy 

0.45 

* 10® 

0 . 35 

Aluminum 

10 

* 10® 

0 . 3 

Steel 

30 

* 10® 

0 . 3 


■'‘The total G values of reference [23] were scaled down by 
factor 2 because Erdogan and Gupta defined G for an isotropic, 
homogeneous plate as 

G - 2 — ^ ( 1 - 2 ) ) 

^Values in parenthesis were obtained with simpler formulas. 
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(a) Normal stress ( Oy) distribution ahead of the crack 


i 


y 



|— A A 

(b) Four 4-noded elements at the crack tip 


Fig. 1: Irwin's virtual crack- closure method 
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o 


y 


(c) Collapsed 4-noded 

elements at crack tip 


(b) Stress distribution and nodal 
forces 


Lg. 2: 


Normal Stress Distribution and Consistent Normal Forces 







(b) QP element idealization (c) Cubic singularity 

element idealization 


Fig. 5 - Forces at and ahead of the crack tip 
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(b) Superposition of computed forces and forces needed 
to erase the pressure loading 


Modification of forces for uniform crack face pressure 
loading 





8-noded parabolic element 


QP singularity element 



12-noded cubic element 


Cubic singular element 


Fig. 7 - Consistent nodal forces for uniform pressure loading in 
various elements co 





(a) CCT specimen 

( h/b » 2.0; a/b *= 0.8 ) 



h a 

(c) Finite element model 
(34 elements) 


Fig. 8 - Mode I problems and 


(b) SEN specimen 

( h/b * 1.0; a/b » 



(d) Detail near the crack tip 
( A/a = 1/16 ) 


finite element idealization 




( h/b « 1.5; a/b « 0.1 
uniform crack face 


(a) Angle crack in a finite plate pressure loading, s) 

(a/b - 0.5) 


Fig. 9 - Mixed mode problems 






Fig. 11: Computation of Normal and Tangential Forces 
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(a) Finite element model 
(140 elements) 



(b) Detail of modeling near 
the crack tip 
A/a = 0,05 


Fig. 12: Finite Element Idealization for the Bimaterial Plate 

Problem 
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(b) Superposition of forces 
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